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Abstract 

We apply set-valued numerical methods to compute an accurate enclosure of the 
rotation number. The described algorithm is supplemented with a method of 
proving the existence of periodic points, which is used to check the rationality 
of the rotation number. A few numerical experiments are presented to show 
that the implementation of interval methods produces a good enclosure of the 
rotation number of a circle map. 

Keywords: Rotation number, Circle map. Rigorous computation, Interval 
arithmetic 


1. Introduction 

Given an orientation-preserving homeomorphism /: 5"^ —> S'^ of the circle, 
the rotation number p is a topological invariant of /, defined as p = p{f) = 
p(F)(mod 1), where 

P(F) := li„ (1) 

n—>-(x> Tl 

Here the function F - a. lift of / - is a homeomorphism of the real line, satisfying 
/ o tt{x) = tt o F{x), where tt: K —>■ 5”^ is the natural projection given by 
Tr{x) = a;(mod 1). By a classical result of Poincare [T] the limit Qis well- 
defined up to an integer, and independent of the choice of x. The circle map / 
has a periodic point if and only if its rotation number is rational: if p = p/g, 
where p and q are co-prime integers, then / has a periodic point of prime period 
q. If, on the other hand, the rotation number is irrational, then there are two 
possibilities: either all orbits are dense in and / is topologically conjugate 
to a rigid rotation or there exists a Cantor set C C which is invariant under 
/ and both forward and backward orbits of all points converge to C, and then 
/ is semi-conjugate to a rigid rotation. 

Apart from a rich flora of theorems concerning rotation numbers for generic 
classes of circle maps, there exists several results concerning the estimation of 
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the rotation number p of a given circle map /, see e.g. m, 0, a- All such 
existing methods require the computation of very long trajectories of /; this 
is apparent from the basic definition Q of the rotation number. When the 
dynamics of / displays sensitive dependence (due to regions of expansion), it is 
well-known that numerical computations of even moderately long trajectories 
can be very inaccurate. As such it may be warranted to question the reliability 
of numerically computed estimates of a rotation number. We address this issue 
by presenting a method that produces explicit, computable bounds on p for a 
given /. As a result, we can e.g. guarantee the correctness of the n first decimals 
of our computed estimate of p. 

In order to guarantee a mathematical rigour in our numerical work, we base 
all computations on set-valued mathematics. This approach is based on interval 
arithmetic [SHIZ], and is very well suited for controlling the propagation of 
numerical errors throughout long iterative processes. 

This paper has the following structure: we begin by recalling some algo¬ 
rithms that - theoretically, and under the right circumstances - can be used 
to produce enclosures of a rotation number. Next, we make the most efficient 
of these algorithms practically applicable by deriving explicit bounds that are 
needed in the computations. The end-product is a tight enclosure of the sought 
rotation number p. By adding a sub-algorithm that proves the existence (and 
uniqueness) of moderately long periodic orbits, we can also handle the case of 
rational rotation numbers. Finally, we present some numerical results from our 
implementation of the described algorithm. 

2. Algorithms for computing the rotation number 

The most simple-minded algorithm for computing the rotation number is a 
direct application of the basic definition 0- 

Theorem 2.1. Let Pn{x) = ^ ■ Then \pn{x) — p\ < for all x and N. 

This method produces explicit bounds of the rotation number; its main draw¬ 
back is that the convergence is only linear. So, given an initial point Xq G S^, 
and an (interval) enclosure a; at of F^{xo), we can define pjy = , which 

produces the bounds: 


P-N ~ <P^Pn + 1 /^- 

Here, we denote the endpoints of an interval according to a; = [x, a:]. 

In [8] the rotation number p is constructed via the continued fraction expan¬ 
sion obtained from the topological behaviour of the circle map: 


P = 


1 


Oi -f 



[0; oi, 02,... ]. 


( 2 ) 
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The continued fraction coefficients at depend on increasingly high iterates of the 
circle map, and - from a numerical point of view - they are not easy to com¬ 
pute. Nevertheless, using this approach, it is possible to construct a sequence of 
intervals of rapidly decreasing lengths, all containing the rotation number p. In 
general the convergence is quadratic, and when the rotation number satisfies a 
Diophantinej^ condition, the convergence can be made cubic, see [H]. Note that 
the set of rotation numbers satisfying a Diophantine condition has full measure 
in [0,1]. On the other hand, the complement lies dense in [0,1]. In what follows, 
we will not assume any Diophantine properties of the rotation number. 

Let / be an orientation preserving homeomorphism on the circle = K/Z. 
Set 6o = /(O) e (0,1), and define /q to be the smallest of the two intervals 
[bo, 1) and (0, 6o]. If bo = 1/2 there is a tie, and we set Iq = (0,1/2]. Next, let 
oi be the first return time of /(&o) to the interval Jq, i.e., the smallest positive 
integer such that f°‘^{bo) G Iq- If there is no such return, we set oi = oo, and 
terminate the process. Otherwise we define 6i = /““^(O), and define Ii to be the 
interval [&i,l) or (0,6i] that is disjoint from Jg. This ends the initial stage of 
the general procedure for computing the continued fraction of p. 

Continuing with the general inductive stage, we let {Ii\ be a sequence of 
half-open intervals of the form [6^,1) or (0, 6i], and let {(fi}, i > I denote the 
sequence of hrst return maps of / from Ii to [Ii IJ li-i). We define to be 
the smallest positive integer such that {bf) G Ii- Given this return time, 

we define As before, if the return map is not well-defined, 

we set Oi+i = 00 and terminate the process. In parallell with this process we 
recursively define positive integers pi and qi as followt]^ 

p_i = 1, Po = 0; p,+i = ai+ip, + pi^i, (3) 

q-i =0, qo = 1; qi+i = ai+iqi + g^-i. 

These integer sequences are used in the following theorem which provides explicit 
bounds on the rotation number p. 

Theorem 2.2. Let Ni be the number of iterates needed to compute qi. For 
the sequences defined above the following holds: 

(a) If p is irrational, then ^ converges to p as i ^ oo. If p is rational, then 
the process terminates (Oi+i = oo) and the last estimate ^ = p. 

(b) If Pi and qi are found, then p is contained in a closed interval A with 

end-points ^ and ; where the integer a is a lower bound of 

O-i+l- 


real number p is Diophantine if there exist positive B and 7 such that \p—p/q\ > 
for all rational numbers p/q. 

^We are assuming that the rotation number is less than 1/2. For p £ [1/2,1) the recursion 
has different initial conditions. 
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(c) \A\ < 4/iVf. For any Ni < N < Ni+i, \A\ < 2/{q^Ni). If 

satisfies the Diophantine condition a^+i < BqJ for some B and 7 > 0, 
then |A| < 2{B + 2 )i/(i+ 7 )(i/JV)^+V(i+-) _ 

We will use a version of case (b). First, let us begin by stating an obvious 
fact about continued fractions of rotation numbers. 

Lemma 2.3. Let oi,..., be the leading elements of the continued fraction of 
the rotation number p: 

p= [ 0 ;ai,a 2 ,...,ai,...] 

then we have the enclosure p € [ 0 ; oi,..., + [ 0 , 1 ]], which is short-hand for the 

interval with lower endpoint [ 0 ; Oi,..., + 1 ] and upper endpoint [ 0 ; Oi,..., a^]. 

Proof. Recall that all are positive integers by construction. The upper 
endpoint [0; Oi,..., a^] corresponds to the case = 00 . The lower endpoint 
is obtained when = 1 and 0^+2 = 00 . 

Remark. If we knew that the rotation number was irrational, we would get a 
tighter bound: p G [0; Ui,..., + [0, cr]], where a = (\/5 — l)/2 « 0.6180339887 

is the golden ratio conjugate. Its continued fraction is cr = [0; 1,1,... ]. 

By using the procedure of Theorem |2.2[ it is clear that we need a means 
to compute very high iterates of the circle map /. This is needed in order to 
compute the return maps (fi, which in turn are used to determine the return 
times Oi appearing in the recursions ([^, and in the continued fraction of p. 
But computing high iterates of a non-linear dynamical system is not straight¬ 
forward. We will address this issue in the following section. 

3. The interval Newton method and shooting techniques 

As already mentioned, the described algorithms all require long pieces of tra¬ 
jectories in order to produce good approximations (and bounds) on the rotation 
number. Naturally, this is problematic as it is not always clear how to control 
the propagation of numerical errors in such computations. Therefore we need a 
reliable and accurate method to compute long trajectories. A common solution 
to this problem is to use multi-precision numerical software. This will allow for 
longer accurate trajectories, but the penalty is high: a decrease in speed by a 
factor 200-1000 is not uncommon. We will illustrate this in Section!^ 

Sticking to standard floating point computations, we will use a version of 
the shooting method described in [9] but modified to suit our particular set-up. 
This method is based on a combination of the interval Newton method and 
multiple shooting. 

The interval Newton method is a constructive implementation of the bounds 
required in Kantorovic’s theorem in order to guarantee the convergence of New¬ 
ton’s method. As such, it can be used to prove the existence (or non-existence) 
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of zeros of general n-dimensional maps. Let g : K" —>■ M" be a continu¬ 
ously differentiable function, and suppose that we have an interval extension 
of its derivative, i.e., given an n-dimensional interval variable jz C K", that is 
2 : = (iCi, ■ • • , Xn), we can compute the interval image Dg{z) 3 {Dg{z): z S z}. 

To verify the existence of zeros of g, we introduce the interval Newton oper¬ 
ator 

N,[z) = z-[Dg{z)]-^g{z), (4) 

where [Dg{z)\ is an interval matrix containing all Jacobian matrices of g of the 
form Dg{z) for z G z. Here z is an arbitrary point from the interval vector 
z usually chosen to be the middle point of z. The operator Q possesses the 
following key properties: 

Theorem 3.1. 

1) If Ng{z) C 2; then there exists exactly one point z* & z such that g{z*) = 0. 

2) If Ng{z) n 2: = 0 then there are no zeros of g in z. 

Let 2:0 = 2; be the initial enclosure of a possible zero of g, and define the 

sequence of intervals 2:^+1 = Ng{zk) n Zk, fc = 0,1,2,_ If a true zero 2* 

is contained in Zq, and if the interval Newton operator is well-defined on this 
domain, then the operator remains well-defined for all iterations, we have 2* G 
Zk, and the intervals Zk form a nested sequence converging to the zero of g. 

By applying the interval Newton operator to T’ - a high-dimensional variant 
of / based on the shooting technique - we can efficiently recast the problem of 
tracking long orbits of / into one about finding zeros of F. To be more precise: 
let f : ^ he a circle map as defined in Section [ 2 ] In order to compute 

the n iterates xi,X 2 , ■ ■ ■ ,Xn, where Xk = f^{xo), we apply the interval Newton 
operator to the n-dimensional map F{-;xo) ■ —>■ (S”^)" defined as follows: 

Fi{z;xo) = f{xo)-xi, 

F2{z;xo) = fixi)-X2 


F„{z;xo) 


f{Xr,-l) - Xn, 


where 2 = (xi, • • • , x„) and xq is the initial point of the orbit. 

Note that the Jacobian matrix of F{-;xo) at 2 = (xi, • • • ,x„) has the fol¬ 
lowing - very sparse - form: 


DF{z]Xo) 


( 1 


0 


0 

1 

-f'{x2) 


0 

0 

1 


0 \ 

0 

0 


V 0 0 


-/'(x„_i) 1/ 


To get a rigorous enclosure of a finite length trajectory starting 
begin by computing an approximate trajectory 2 = (xi,-- - ,x„). 


at xo, we 
We then 
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consider a neighbourhood of this piece of trajectory by widening each component 
into an interval, producing the box z = {xi, ■ ■ ■ , Xn). Next, if we can show that 
the interval Newton operator 

Nf{z] xo) = z- DF{z; Xo)~^F{z-, Xq) 

maps this interval vector into itself Nf{z;xo) C z, then we have a true trajec¬ 
tory enclosed in 2 :. By reapplying the interval Newton operator a few times, we 
can tighten the enclosure of the trajectory as much as we want. This gives us 
very high-quality trajectories, starting from only coarse approximate trajecto¬ 
ries. 

The computation of the interval matrix inverse can cause difficulties for large 
n. Fortunately, the explicit inverse is not needed; we can simply solve for the 
correction term h = DF{z-,xo)~^F{z;xo). This equation can be recast as 

DF{z;xo)h = F{z;xo), 


where z = {xi, ■ ■ ■ ,Xn), z = (xi, • • • , x„) and h = {hi, • • • , hn). 

We define all i recursively in the following way 

hi = f{xo)-xi, 

h2 = f{xi)-X2 + f{xi)hi, 

hn — f{Xn—l) Xji f {Xn—l)hn—l~ 

On completion of these computation, we can evaluate the interval Newton op¬ 
erator as Nf{z; Xq) = z — h. This concludes the description of how to compute 
long, yet highly accurate (in fact rigorous), iterates of the circle map /. 


4. Verification of the rationality of the rotation number 


Recall from part (a) of Theorem 2.2 that if p is rational, then the process 


terminates (ci+i = 00 ) and the last estimate ^ = p. In fact, this means that 
the circle map / has a stable period-gi point, and all forward trajectories tend 
to the corresponding periodic orbilQ With this in mind, we always search for a 
period-gi point of the circle map before attempting to determine the next return 
time Oi+i. 

Again we use the interval Newton method to verify the existence of a periodic 
point. The line of argument is very similar to that of the previous section. We 
first try to locate a good candidate by iterating a random point many times, 
producing the finite orbit xq,Xi, ... ,Xn■ Here N is taken much larger than the 
sought period g^. If there is a period orbit, these iterates should accumulate on 
it, which means that xn+i,xn+ 2 ^ ■ ■ ■ ^XN+qi should form a good approximation 


®The periodic point is attracting. 
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of the periodic orbit. Widening this finite set of points into a gi-dimensional 
box, and verifying that this box is mapped into itself under the interval Newton 
map, established the existence of a period-gi point. This - in turn - means that 
ai_|_i = oo, and that the rotation number is rational: p = ^. 


5. Numerical experiments 

In this section we present numerical results from our implementations of both 
described algorithms: the linear one (based on Theorem 2.1) and the quadratic 
one (based on Theorem |2.2[ ). We base our code on the CAPD interval library 
|10j . For all computations performed, we list the timings in seconds using a 
single thread on an Intel(R) Core(TM) i7 CPU 920 @ 2.67GHz. 


5.1. The Arnold family 

We consider the two-parameter Arnold family of circle maps given by 
fa,e{x) = X + a — esin(27ra;)(mod 1). 

The map fa,e{x) with 27r|e| < 1 is an orientation-preserving analytic circle 
diffeomorphism for every pair of parameters (a, e). The set of parameters (a, e) 
corresponding to a fixed irrational p forms a continuous curve known to be 
analytic when p is a Diophantine number; for the rational p the set has an non¬ 
empty interior, and formes the well-known the Arnold tongue of the rotation 
number p El [Ill- 

Table [2 displays enclosures of the rotation number when computed according 
to Theorem |2.1| using set-valued computations for different pairs of parameters 
of the Arnold map. We list the final enclosure of the rotation number, its 
radius, and the number of iterations of the map /. In order to save space, and 
to highlight the significant digits, we use a compact notation for intervals: for 
example, 0.2fg means the interval [0.218,0.221]. As expected, our numerical 
results indicate that this method gives only linear convergence. 


N 

P 1 rad(p) 

P 

rad(p) 


a = 0.22 e = 0.01 

a = 0.45 e = 0.01 

1000 

5000 

10000 

15000 

0.2^^ 

0.2?° 

0.219f} 

0.219fl 

1.0 X 10"^ 
2.0 X 10-4 
1.0 X 10-4 
6.7 X 10-9 

^^9 

o.4t^i 

1.0 X 10-9 
2.0 X 10-4 
1.0 X 10-4 
6.7 X 10-9 


a = 0.22 e = 0.159 

a = 0.45 e = 0.159 

1000 

5000 

10000 

15000 

Iklfifl 

0.162^1 

0.162^ 

0.16249 

1.0 X 10-9 
2.0 X 10-4 
1.0 X 10-4 
6.7 X 10-9 

^(461 

0.4649 

0.461448 

0.46149 

1.0 X 10-9 
2.0 X 10-4 
1.0 X 10-4 
6.7 X 10-9 


Table 1: Results for the Arnold map; linear convergence. 


7 



















N 

prec (bits) 

T (sec) 

prec (bits) 

T (sec) 


a = 0.22 

e = 0.01 

a = 0.45 

e = 0.01 

1000 

1000 

0.06 

1000 

0.05 

5000 

1000 

0.27 

1000 

0.29 

10000 

1000 

0.56 

1000 

0.56 

15000 

1000 

0.84 

1000 

0.83 


a = 0.22 

e = 0.159 

a = 0.45 

e = 0.159 

1000 

9500 

1.83 

10500 

2.25 

5000 

9500 

9.20 

10500 

11.24 

10000 

9500 

18.42 

10500 

22.34 

15000 

9500 

27.82 

10500 

33.82 


Table 2: Timings in seconds for the algorithm with linear convergence for the 
Arnold map. 


Enclosures of the rotation number obtained using our algorithm (Lemma |2.3| ) 
are given in Table Here the convergence is at least quadratic, and the use of 
interval arithmetic allows us to get a much better enclosure already for a small 
number of iterations. 



P 

rad(/o) 


P 

rad(p) 

i 

a = 0.22 

e = 0.01 

i 

a = 0.45 

e = 0.01 

4 

r\ 

U.^195 

2.4 X 10"'‘ 

3 

O 

OOC 

8.6 X 10’^ 

6 

0.219f| 

3.9 X 10"® 

4 

0.449975?? 

1.2 X 10“'^ 

8 

0.2198J? 

1.6 X 10"® 

4 

0.449975it 

6.1 X 10“® 

9 

0.219809993lf? 

1.1 X 10"“ 

5 

0.4499753532^? 

3.4 X 10““ 

i 

a = 0.22 

e = 0.159 

i 

a = 0.45 

e = 0.159 

3 

“ole^E 

5.0 X 10"® 

3 

“OAi 

1.6 X 10“® 

5 

0.1624iJ 

8.1 X 10"® 

4 

0.461(,J 

1.7 X 10“® 

6 

0.16242i| 

1.4 X 10"® 

5 

0.4610811 

8.7 X 10“'^ 

5 

0.16242it 

5.8 X 10"'^ 

6 

0.4610864?^ 

7.8 X 10“^ 


Table 3: Results for the Arnold map; quadratic convergence. Here i denotes the 
number of coefficients of the continued fraction expansion. 


Using the interval Newton method as described in Section we can often 
verify if the rotation number is rational. For the cases shown in Table no 
periodic points were detected, and so the results encloses a (possibly) irrational 
rotation number. Similar computations were performed for other pairs of pa¬ 
rameters (a, e) given in Table These parameters produce rational rotation 
numbers, and thus we can provide an exact representation of the form p = |, 
where p and q correspond to the integers pi, qi in Theorem 2.2 Note that here 
q is the period of the periodic point. 

We end this example by computing rotation numbers for many different pa¬ 
rameters. This produces the well-known devil’s staircase graph, illustrating the 
frequency locking phenomena. To be concrete, we will take e G {0.01, 0.1,0.159}, 



























a 

€ 

T (sec) 

0.22 

0.01 

2.99 

0.22 

0.159 

2.57 

0.45 

0.01 

40.88 

0.45 

0.159 

1.82 


Table 4: Timings in seconds for the algorithm with quadratic convergence for the 
Arnold map. The large timing for (a, e) = (0.45, 0.01) is due to a large coefficient 
(04 = 100 ) appearing in the continued fraction of the rotation number. 


a 

e 

P 

p/q 

e) 

T (sec) 

0.22 

0.159154943 

0.162 

6/37 

0.013173605 

0.62 

0.2 

0.159 

0.125683060 

23/183 

0.012343856 

0.95 

0.21 

0.159 

0.142857143 

1/7 

0.037952001 

0.61 

0.43 

0.159154943 

0.428571 

3/7 

0.021337727 

0.61 

0.21 

0.15 

0.153846154 

2/13 

0.005169938 

0.23 

0.16 

0.1589 

0.018957346 

4/211 

0.856340865 

0.78 


Table 5: Rational rotation numbers for the Arnold map, and their corresponding 
periodic points x*(a,e). 


and for each value of e we compute the enclosure of the rotation number p for 
1000 different values of a. The results are illustrated in Figure 



Figure 1: The devil’s staircase for e £ {0.01,0.1,0.159}. 
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Figure 2: Successive zooms of the devil’s staircase for e = 0.159. 


5.2. The delayed logistic map 

Similar computations were carried out for the delayed logistic map, defined 

The approximation of the rotation number for this map has been considered 
in [T^. It is believed that for 2 < A < 2.16 the map (j) possesses a smooth 
invariant curve. In order to prove the existence of an invariant curve for the 
delayed logistic map, one may adopt the method described in m. where a 
topological approach to prove the existence of a normally hyperbolic manifold 
for maps is introduced. As this is not the focus of our work, we will simply 
assume that the invariant curve exists for the parameters we are considering. 


A 

Pv 

P 

rad(p) 

T (sec) 

2.04 

0.1628037 


2.4 X 10"^ 

0.003 

2.06 

0.1607109 

0.160111 

1.1 X lO"-* 

0.004 

2.08 

0.1584864 

0.158|| 

6.0 X 10-5 

0.005 

2.10 

0.1561058 

0.15ii 

3.6 X 10-* 

0.005 

2.12 

0.1535363 

0.15351? 

8.0 X 10-6 

0.005 

2.14 

0.1507185 

0.150i| 

1.3 X 10-* 

0.005 

2.16 

0.1474935 

0.147|? 

8.6 X 10-5 

0.005 


Table 6 : Computational results for the delayed logistic map, using the algorithm 
with quadratic convergence. The short computational times are due to the 
simple form of the map (no trigonometric functions). 

To compute an approximate enclosure of the rotation number, we first iterate 
the delayed logistic map several times in order to avoid transient behaviour. 
This leaves us with a point very close to the invariant curve (see Figure [^, 
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and we will use this point as our initial point. To compute intervals \Ii\ from 
Theorem 2.2 at each iteration we consider an angle defined using arctangent 
function of two arguments. We present our results for different values of the 
parameter A in Table For comparison, the second column py was computed 
using the method presented in |12j . 



Figure 3: The invariant curves for the delayed logistic map for A € {2.05, 2.10, 2.15}. 


6. Discussion 

We have presented a method for computing an enclosure of the rotation 
number of a given circle map. This method is based upon set-valued compu¬ 
tations combined with explicit bounds that can be obtained whilst computing 
the continued fraction of the sought rotation number. We have demonstrated 
our method by applying it to the classical two-parameter Arnold family of circle 
maps, and to the delayed logistic equation. In the latter case, we compare our 
results with those appearing in [T^ . 

Of course, our main goal so far has been to obtain mathematical rigour in 
the numerical computations necessary. Nevertheless, by utilizing the method 
of multiple shooting together with the interval Newton method, we can per¬ 
form all computation in double precision. This makes our method quite fast in 
comparison to approximative methods that utilize multi-precision. 
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